Setup

Install any packages you might be missing:

library(spatstat)
library(splancs)
# tree data
data(bramblecanes)
data(finpines)
data(spruces)

#get info
print(finpines)
## Marked planar point pattern: 126 points
## Mark variables: diameter, height 
## window: rectangle = [-5, 5] x [-8, 2] metres
intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern:  126 points
## Average intensity 1.26 points per square metre
## 
## Coordinates are given to 13 decimal places
## 
## Mark variables: diameter, height
## Summary:
##     diameter         height     
##  Min.   :0.000   Min.   :0.800  
##  1st Qu.:1.000   1st Qu.:1.825  
##  Median :2.000   Median :2.850  
##  Mean   :2.532   Mean   :2.828  
##  3rd Qu.:3.000   3rd Qu.:3.600  
##  Max.   :7.000   Max.   :5.400  
## 
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
#creates point pattern data format from the x and y coordinates
finpines.pts<-as.points(finpines$x,finpines$y)
bramblecanes.pts <-as.points(bramblecanes$x,bramblecanes$y)
spruces.pts <-as.points(spruces$x,spruces$y)

Exploring patterns

dev.off()
## null device 
##           1
#divide into quadrants and count number of points in the quadrant
Q <- quadratcount(finpines, nx = 4, ny = 3)
Q
##                x
## y               [-5,-2.5) [-2.5,0) [0,2.5) [2.5,5]
##   [-1.33,2]            10        8      11      21
##   [-4.67,-1.33)         8        7       9       3
##   [-8,-4.67)           12       11      22       4
plot(finpines.pts,pch=19, main="Finpines", xlab="Easting",ylab="Northing")
plot(Q, add = TRUE, cex = 2)
# Finpines data
sp_point <- matrix(NA, nrow=length(finpines$x),ncol=2)
sp_point[,1] <- finpines$x
sp_point[,2] <- finpines$y
colnames(sp_point) <- c("x","y")
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="Easting",ylab="Northing",cex=.5)

intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern:  126 points
## Average intensity 1.26 points per square metre
## 
## Coordinates are given to 13 decimal places
## 
## Mark variables: diameter, height
## Summary:
##     diameter         height     
##  Min.   :0.000   Min.   :0.800  
##  1st Qu.:1.000   1st Qu.:1.825  
##  Median :2.000   Median :2.850  
##  Mean   :2.532   Mean   :2.828  
##  3rd Qu.:3.000   3rd Qu.:3.600  
##  Max.   :7.000   Max.   :5.400  
## 
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
# Random points
u.x <- runif(n=nrow(sp_point), min=bbox(sp_point)[1,1], max=bbox(sp_point)[1,2])
u.y <- runif(n=nrow(sp_point), min=bbox(sp_point)[2,1], max=bbox(sp_point)[2,2])

# Regular points
r.x <- seq(from=min(sp_point[,1]),to=max(sp_point[,1]),length=sqrt(nrow(sp_point)))
r.y <- seq(from=min(sp_point[,2]),to=max(sp_point[,2]),length=sqrt(nrow(sp_point)))
r.x <- jitter(rep(r.x,length(r.x)),.001)
r.y <- jitter(rep(r.y,each=length(r.y)),.001)

# Plot Finpines, uniform, and regular points
par(mfrow=c(1,3),mar=c(4,4,1.5,0.5))
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="Easting",ylab="Northing",cex=.5,pch=19)
plot(x=u.x,y=u.y,main="Random Points", xlab="Easting",ylab="Northing",cex=.5,pch=19)
plot(x=r.x,y=r.y,main="Regular Points", xlab="Easting",ylab="Northing",cex=.5,pch=19)

Intro to homogeneous Poisson processes

# CSR from random uniform point process
## Generates a CSR point process conditionally on a fixed number N of points

## Let N = 100
par(mfrow=c(2,2),mar=c(0.4,0.4,1,0.4))
pp.100.unif<-runifpoint(100,win=owin(c(0,1),c(0,1)))
pp.100.unif
## Planar point pattern: 100 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unif)
## [1] 100
plot(pp.100.unif, main="")
title('100 uniform points')


# repeat a few times to see different realizations of the process

# HPP with intensity 100 in the unit square
pp.100.unit <- rpoispp(100,win=owin(c(0,1),c(0,1)))
pp.100.unit
## Planar point pattern: 107 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unit) 
## [1] 107
plot(pp.100.unit,main="")
title('HPP intensity = 100, unit square')

# repeat a few times to see different realizations of the process (including different numbers of points)

# HPP with intensity 1 in a 10 x 10 square
pp.1.10 <- rpoispp(1, win=owin(c(0,10),c(0,10)))

par(mfrow=c(1,2))
plot(pp.100.unit,main="")
title('intensity = 100, unit square')
plot(pp.1.10,main="")
title('intensity = 1, 10 x 10 square')

# HPP with intensity 200 in the unit square
par(mfrow=c(1,1))
pp.200.unit <- rpoispp(200)
plot(pp.200.unit)

# HPP with intensity 100 in the unit square
side.square <-1
pp.100.unit <- rpoispp(10000,win=owin(c(0,side.square),c(0,side.square)))
points.tot <- npoints(pp.100.unit) 
plot(pp.100.unit,main="")

# inscribe a circle within the square
radius.circle <-  side.square/2
win.circle <- disc(radius.circle, c(radius.circle,radius.circle))

pp.incircle <- pp.100.unit[win.circle]
points.incircle <- npoints(pp.incircle)
summary(pp.incircle)
## Planar point pattern:  7942 points
## Average intensity 10116.13 points per square unit
## 
## Coordinates are given to 8 decimal places
## 
## Window: polygonal boundary
## single connected closed polygon with 128 vertices
## enclosing rectangle: [0, 1] x [0, 1] units
## Window area = 0.785083 square units
## Fraction of frame area: 0.785
plot(pp.incircle)

# calculate pi using ratio of areas and number of points
piEst <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
piEst
## [1] 3.13913

Function to simulate pi over many iterations

pi.simulation2 = function(iter, lambda, side.square){
  
  piEst.list = rep(0,iter)
  points.incircle.list = rep(0,iter)
  
  # inscribe a circle within the square
  radius.circle <-  side.square/2
  win.circle <- disc(radius.circle, c(radius.circle,radius.circle))
  
  for(i in 1:iter){
    # HPP with intensity 100 in the unit square
    pp.square <- rpoispp(lambda,win=owin(c(0,side.square),c(0,side.square)))
    points.tot <- npoints(pp.square) 
    
    # find points within the circle
    pp.incircle <- pp.square[win.circle]
    points.incircle <- npoints(pp.incircle)

    # calculate pi using ratio of areas and number of points
    pihat <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
    piEst.list[i] = pihat
  }
  return( piEst.list)
}

iter = 1000
lambda = 1000
side.square = .5
est = pi.simulation2(iter,lambda, side.square)
mean(est)
## [1] 3.144389
var(est)
## [1] 0.01097336
# we plot the histogram and see it appears to be binomial. why?
hist(est,nclass = 15)

Tests for CSR

# create polygon (square) area for finpines data
minfx <- min(finpines$x)
maxfx <- max(finpines$x)
minfy <- min(finpines$y)
maxfy <- max(finpines$y)

polyfin <- as.points(c(minfx,maxfx,maxfx,minfx),c(minfy,minfy,maxfy,maxfy))
UL.khat.fin <- Kenv.csr(length(finpines$x), polyfin, nsim=99, seq(0.05,3.5,0.15))
## Doing simulation  1 
## Doing simulation  2 
## Doing simulation  3 
## Doing simulation  4 
## Doing simulation  5 
## Doing simulation  6 
## Doing simulation  7 
## Doing simulation  8 
## Doing simulation  9 
## Doing simulation  10 
## Doing simulation  11 
## Doing simulation  12 
## Doing simulation  13 
## Doing simulation  14 
## Doing simulation  15 
## Doing simulation  16 
## Doing simulation  17 
## Doing simulation  18 
## Doing simulation  19 
## Doing simulation  20 
## Doing simulation  21 
## Doing simulation  22 
## Doing simulation  23 
## Doing simulation  24 
## Doing simulation  25 
## Doing simulation  26 
## Doing simulation  27 
## Doing simulation  28 
## Doing simulation  29 
## Doing simulation  30 
## Doing simulation  31 
## Doing simulation  32 
## Doing simulation  33 
## Doing simulation  34 
## Doing simulation  35 
## Doing simulation  36 
## Doing simulation  37 
## Doing simulation  38 
## Doing simulation  39 
## Doing simulation  40 
## Doing simulation  41 
## Doing simulation  42 
## Doing simulation  43 
## Doing simulation  44 
## Doing simulation  45 
## Doing simulation  46 
## Doing simulation  47 
## Doing simulation  48 
## Doing simulation  49 
## Doing simulation  50 
## Doing simulation  51 
## Doing simulation  52 
## Doing simulation  53 
## Doing simulation  54 
## Doing simulation  55 
## Doing simulation  56 
## Doing simulation  57 
## Doing simulation  58 
## Doing simulation  59 
## Doing simulation  60 
## Doing simulation  61 
## Doing simulation  62 
## Doing simulation  63 
## Doing simulation  64 
## Doing simulation  65 
## Doing simulation  66 
## Doing simulation  67 
## Doing simulation  68 
## Doing simulation  69 
## Doing simulation  70 
## Doing simulation  71 
## Doing simulation  72 
## Doing simulation  73 
## Doing simulation  74 
## Doing simulation  75 
## Doing simulation  76 
## Doing simulation  77 
## Doing simulation  78 
## Doing simulation  79 
## Doing simulation  80 
## Doing simulation  81 
## Doing simulation  82 
## Doing simulation  83 
## Doing simulation  84 
## Doing simulation  85 
## Doing simulation  86 
## Doing simulation  87 
## Doing simulation  88 
## Doing simulation  89 
## Doing simulation  90 
## Doing simulation  91 
## Doing simulation  92 
## Doing simulation  93 
## Doing simulation  94 
## Doing simulation  95 
## Doing simulation  96 
## Doing simulation  97 
## Doing simulation  98 
## Doing simulation  99
khat.fin<-khat(finpines.pts, polyfin, seq(0.05,3.5,0.15))
# plot of Khat-pi t^2 from data
plot(seq(0.05,3.5,0.15), khat.fin-pi*seq(0.05,3.5,0.15)^2, type="l", xlab="Distance", ylab="Estimated K-pi h^2")

# plot upper bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$upper-pi*seq(0.05,3.5,0.15)^2, lty=2)

# plot lower bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$lower-pi*seq(0.05,3.5,0.15)^2, lty=2)

# L functions
l<-function(k,h)
{sqrt(k/pi)-h}

# plot Lhat from data
plot(seq(0.05,3.5,0.15), l(khat.fin, seq(0.05,3.5,0.15)), type="l", xlab="Distance", ylab="Estimated L",ylim=c(-.1,.2))

# plot upper bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$upper,seq(0.05,3.5,0.15)), lty=2)

# plot lower bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$lower,seq(0.05,3.5,0.15)), lty=2)

K<-Kest(finpines)
plot(K)

plot(K, cbind(r, sqrt(iso/pi)) ~ r)

plot(K, cbind(trans,iso,border) - theo ~ r)

plot(envelope(finpines,Kest))
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

# G test: Nearest neighbor distribution test
Gtest<-Gest(finpines)
plot(envelope(finpines,Gest))
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

# F test: Point to nearest event distribution test
Ftest<-Fest(finpines)
plot(envelope(finpines,Fest))
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

# J test: J(r) = (1-G(r)/(1-F(r)))
Jtest<-Jest(finpines)
plot(envelope(finpines,Jest))
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

Inhomogeneous Poisson Process

lamb1 <- function(x,y) {100 * exp(10*x-5*y)}
Ipp1 <- rpoispp(lamb1, 100)
Ipp1.p<- as.points(Ipp1$x,Ipp1$y)


plot(Ipp1.p,main="IPP, intensity(x,y)=100*exp(10x-5y)",pch=20,xlab='x',ylab='y')

lamb2 <- function(x,y) {100 * exp(-10*x+5*y)}
Ipp2 <- rpoispp(lamb2, 100)
Ipp2.p<- as.points(Ipp2$x,Ipp2$y)


plot(Ipp2.p,main="IPP, intensity(x,y)=100*exp(-10x+5y)",pch=20,xlab='x',ylab='y');

# can save a step and write function within rpoispp
Ipp3 <- rpoispp(function(x,y) {100 * exp(-3*x)}, 100)
plot(Ipp3,pch=20,xlab='x',ylab='y')

  • Simulations of (2) clustered point processes
#Generate a square polygon
spoly<-as.points(c(0,1,1,0),c(0,0,1,1)) #polygon

#intensity of parents is 15, average number of offspring is 7, and variance of distance traveled by child (away from parent) in 0.0025
pcp1<-pcp.sim(15,7,0.00025,spoly)
plot(pcp1,main="PCP, (P,O,Spread)=(25,7,.00025)",pch=20,xlab='x',ylab='y')

#intensity of parents is 25, average number of offspring is 4, and variance of distance traveled by child in 0.005
pcp2 <-pcp.sim(25,4,0.005,spoly)
plot(pcp2,main="PCP, (P,O,Spread)=(25,4,.005)",pch=20,xlab='x',ylab='y')

Estimating the parameters of a cluster process

# Using K-function for Clustered Process
# Simulate homogeneous pcp (same as above)
pcp2 <-pcp.sim(25,4,0.005,spoly)
# need to define polygonal area
spoly<-as.points(c(0,1,1,0),c(0,0,1,1))
# estimate
pcp.fit<-pcp(pcp2,spoly,h0=.3,n.int=40)

#simulate from the theoretical process with parameter estimates from above

plot(pcp2)
m.h <- npts(pcp2)/(areapl(spoly)*pcp.fit$par[2])
sims.h <- pcp.sim(pcp.fit$par[2], m.h, pcp.fit$par[1], spoly)
pointmap(as.points(sims.h), add=TRUE, col="red")

# L-function
l<-function(k,s)
{sqrt(k/pi)-s}

# K function for PCP
r <- seq(0.01,1,by=.2)
K.env <- Kenv.pcp(pcp.fit$par[2], m.h, pcp.fit$par[1], spoly, nsim=20, r=r)
L.env<- lapply(K.env,FUN=l,r)

limits <- range(unlist(L.env))
plot(r, sqrt(khat(as.points(pcp2),spoly,r)/pi)-r, ylim=limits, main="L function with simulation envelopes and average", type="l", xlab="distance", ylab="")
lines(r, L.env$lower, lty=5)
lines(r, L.env$upper, lty=5)
lines(r, L.env$ave, lty=6)
abline(h=0)

Application to a real data set (Cardiff)

data(cardiff)

cardiff.pts <-as.points(cardiff$x,cardiff$y)


polymap(cardiff$poly)
pointmap(as.points(cardiff), add=TRUE)
title("Locations of homes of 168 juvenile offenders")

pcp.fit.cardiff <- pcp(as.points(cardiff), cardiff$poly, h0=20, n.int=30)
pcp.fit.cardiff
## $par
##         s2        rho 
## 6.32912666 0.01153204 
## 
## $value
## [1] 0.02499235
## 
## $counts
## function gradient 
##       77       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
m <- npts(as.points(cardiff))/(areapl(cardiff$poly)*pcp.fit.cardiff$par[2])
r <- seq(2,30,by=2)
K.env <- Kenv.pcp(pcp.fit.cardiff$par[2], m, pcp.fit.cardiff$par[1], cardiff$poly,
                  nsim=20, r=r)
L.env<- lapply(K.env,FUN=l,r)

limits <- range(unlist(L.env))
plot(r, sqrt(khat(as.points(cardiff),cardiff$poly,r)/pi)-r, ylim=limits,
     main="L function with simulation envelopes and average", type="l",
     xlab="distance", ylab="")
lines(r, L.env$lower, lty=5)
lines(r, L.env$upper, lty=5)
lines(r, L.env$ave, lty=6)
abline(h=0)